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Abstract. We present a self-consistent analytical model for the computation of the physical conditions in a steady 
quasi-Keplerian accretion disk. The method, based on the thin disk approximation, considers the disk as concentric 
cy linders in which we treat the vertical transfer as in a plane-parallel medium. The formalism generalizes a work 
bv lHubenv I l|l990h . linking the disk temperature distribution to the local energy dissipation and leads to analytical 
formulae for the temperature distribution which help to understand the behaviour of the radiation propagated 
inside the disks. One of the main features of our new model is that it can take into account many heating sources. 
We apply the method first to two sources: viscous dissipation and stellar irradiation. We show that other heating 
sources like horizontal transfer or irradiation from the ambiant medium can also be taken into account. Using the 
analytical formulation in the case of a modified Shakura & Sunyaev radial distribution that allow the accretion 
rate to be partly self-similar in the inner region, and, for an a and l3 prescription of the viscosity, we obtain 
two-dimensional maps of the temperature, pressure and density in the close environment of low mass young 
stars. We use these maps to derive the observational properties of the disks such as spectral energy distributions, 
high resolution spatial images or visibilities in order to underline their different behaviours under different input 
models. 
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1. Introduction 

Since the initial models o f vis cous accretion disks of 
IShaknra, fc Sunvaev I JiflT.Sll and iFvnden-Rel] Pringle I 
l|l974l) . many efforts have focused on the description of 
the close environment of T Tauri stars (hereafter TTS). 
In order to reproduce high angular resolution images 
that became available soon after 1990, various authors 
have published model s of scattered light in large cir- 
cumstellar structures iLazareff. Monin. fc Pudritd fl990t 



IWhitnev & Ha,rtma,nnlll992l) . In the last five years, more 
precise images of edge-on disks have been obtained by 
the Hubble Space Telescope (HS T) and infrared adap- 
tive optics JBurrows et all 11994 [gtapelfeldt et al. 1998; 
iMonin fc Bouvieill200r ft. allowing us to refine the models 
and the physical parameters used to describe the disks. 
Most of the time, when interpreting these images, the au- 
thors use ad-hoc power laws to extract the physical con- 
ditions in the disks, like surface density or the disk scale 
height. Indeed, these images show that the large-scale fiar- 
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ing disks have a complex vertical structure that must be 
taken into account in the interpretation of their flat SED. 

Up to now, we have access to only a very few con- 
straints on the very central parts of circumstellar disks, 
since current instruments do not allow us to directly ob- 
serve within 10 AU of the central object. However, the 
physics of the inner part of the disks is not expected to be 
completely different from that in their outer parts, but our 
knowledge of these central regions of circumstellar disks 
is restricted to scarce interferometric data and indirect 
measurements of the inner phenomena via the spectral 
energy distribution (hereafter SED), the high energy UV 
and X-ray emission of disks or the measurement of mag- 
netic flelds. 

The first models used quasi-Keplerian steady accre- 
tion disks, expected to be geometrically thin for reason- 
able values of the accretion rate. They describe the disk as 
infinitely fiat, predicting that the emergent fiux has a spec- 
tral energy distribution proportional to the power -4/3 of 
the wa velength in the infrared. However. iR.vdgren fc Zaki 
showed that most TTS present fiatter spectral en- 
ergy distributions, i.e. decreasing less steeply or even re- 
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maining constant ov er a large domain in the in frared. 
This discrepancy led I Adams. Shu fc Ladal l|l988l) to as- 
sume that the radial temperature distribution in T Tauri 
disks is also flatter than in standard accretion disks: if one 
assumes that the disk thickness flares then the disk inter- 
cepts more light from the star at large distances and leads 
to a variety of SED slopes. 

Modeling of the vertical structure of the disks 
and the associated r adiati ve transfer has been devel- 
oped by focusing their attention on 
the vertical structure of disks in FU Orionis systems 
where they explained a possibl e origin of outbursts. 
IChiang fc Goldreichl Jl 997l ll flflflh used a semi-analytical 
model of the circumstellar structure to e xplain the SEDs 
of different t ypes of young stars. Recently. IP'Alessio et alJ 
numerically resolved a complete set of equa- 
tions to compute the spatial distribution of the tempera- 
ture and the density that can be constrained by the ob- 
served SEDs and two-dimensional images. 

Up to now, the models have included various approxi- 
mations : (i) the magnetic field is ignored; (ii) the viscos- 
ity is assumed to depend only on local conditions and (iii) 
the viscosity is expressed using a simple ad-hoc law. On 
the other h and, when authors take into account the mag- 
netic field JShn et a,]. IJl Qfl4 iFerreira, fc Pelletier I Fflfli 
ICasse & Ferreiralbnnnalbfl. thev assume that these mag- 
netic fields do no t affect the ove r all dis k str ucture. As 
another example, iTerauera et al. I and iTera uem I 

showed that gravitational instabilities in the disks 
may produce waves dissipating viscous energy far from 
where they are created. Integrating all these effects in a 
unique model appears to be a difficult task. 

Our goal is to develop a radiative model simple enough 
that it could later on include as many detailed physical 
processes of the energy transport in the disks as possible, 
and / or the interaction between the circumstellar mat- 
ter and the magnetic field. In this paper, we present a 
new method to process the radiative transfer i n circum- 
stellar disks, following the pioneering work by iHubenvl 
il99nl) and an initial applica t ion in the case of TTS disks 
hvlMalbet fc RertoutHl99ll):lMa]betlirm9i. Our results 
are not fundamentally different from previous models but 
thanks to the analytical expression of the vertical temper- 
ature distribution that we develop, we are able to iden- 
tify the actual origin of the ten rperature rise in the uppe r 
part of the disk photosphere ifMalbet fc Bertout"lll99ll) . 
We propose that this formaHsm can be the starting point 
for more elaborate models, including magnetic field or vis- 
cous transport. As a first step, our new self-consistent 
analytical model can be used to constrain the physical 
environment of young stars using images from EST or 
ground-based large telescopes, as well as miUimeter radio- 
interferometer maps and forthcoming data from large op- 
tical interferometers , and especially the ESO VLTI project 
ijMalbet et alJl2nnnll . 

Section |2l presents the analytical model derived from 
hydrostatic and radiative equihbrium, as first described 
bv lHubeiivl I^199fl.l . and extended here to the illumination 



by an external source. Sectional describes the application 
of this model to the specific case of the circumstellar en- 
vironment of young stellar objects, especially in the case 
where the accretion rate is not constant at all radii. To 
illustrate the model and its application, we present and 
discuss in Sect. 01 some results in the case where the heat- 
ing by the central source can be neglected and we focus 
our attention on the production of observable quantities 
like SEDs, images, and interferometric data including vis- 
ibility amplitudes and closure phases. 

2. Analytical model 

In this section, we detail the analytical model for the calcu- 
lation of the spatial distribution of quantities like the tem- 
peratu re or the density. We use the notations of lHubenvl 
lll99nh as much as possible, even if we present the formal- 
ism in a slightly different way. The first 3 subsections 112. 11 
12.21 and l2.3|l are short reminders of the main equations of 
the problem. The last one Il2.4|l is the generalization of 
this formaHsm to external sources of radiation. 

The principle for computing the vertical structure of 
T Tauri disks is to consider the disk as a set of concentric 
cylinders of infinitesimal width. Each cylinder is itself re- 
garded as a plane-parallel atmosphere of mixed dust and 
gas layers. The basic assumption is that the radial radia- 
tive transfer of energy is smaller than the vertical transfer 
which is true only in the thin disk approximation l|Pringle I 
However, we can still take into account the horizon- 
tal transfer as a perturbation for disks of moderate fiaring 
(see in Sect. 12.4. Ijl . We solve the coupled equations of hy- 
drostatic equilibrium and the radiative transfer in each 
cylinder in order to compute the temperature and den- 
sity distributions along the z-axis at a given radius r. In 
the following sections, when we focus our attention on the 
vertical structure of the cylinder of radius r, we drop the 
r coordinate unless otherwise specified. We also use the 
mass column coordinate m in place of the height coordi- 
nate z: 



m{z) 



p{z)dz, 



(1) 



where M = to(0) = S/2 is half the surface density of the 
disk at radius r. 



2.1. Initialization stage from the radial structure 

For the initial step, we assume a standard geometrically 
thin disk with a uniform vertical structure at each radius. 
To fully solve the radial structure, we need: 

- the mass per area unit S(r') 

- the energy dissipation per area unit F{r) 

- the opacity law x(T, p) of the disk material 

- the equation of state of the material, p = S{P, T) 

- the local vertical acceleration field Qz (r, z) , generally 
inferred from the gravity 
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where the temperature, density, scale height, pressure and 
optical thickness laws, respectively Td, ^b, Pu, hu and td, 
are given by the coupled equations: 

PDW--g.(r,/iD(r))l](r) 

hB{r) = Y.{r)/pD{r) 
TD(r) = E(r) x(7b(r), Pd W) 



(2) 



In the optically thin case, an iterative computation is 
required. One should notice that the mass, energy dissi- 
pation and acceleration field distributions do not need to 
be independent of the temperature or density laws; in this 
case ^{r), F{r) and Qz{r,z) are refined at each iteration. 

A first refinement is to replace the uniform distribu- 
tion jOD(r) by the isothermal density distribution. The 
resulting radial structure is t hen the structure given by 
the so-called standard model l|Shakura fc Sunvaev lllOT^ 
iLvnden-Bell fc Prinde lll974 . 

Another improvement consists of replacing the effec- 
tive temperature by the temperature of the equatorial 
layer located at rD/2. As in a stellar atmosphere, we sub- 
stitute for the first equation of the set of equations ^ the 
following equation: 



3/4 F(r) (^(r) +2/3) 

1 - c-^d('~) 



(3) 



Theoretically speaking, the radiative transfer equa- 
tions can be applied to any kind of radiation. However, 
common and well-known approximations make them eas- 
ier to solve provided the radiation is thermalized or 
isotropic. For instance we can approximate the fiux- 
weighted mean opacity by the Rosseland opacity if the 
radiation is thermalized with the medium. The Eddington 
closure relations apply well for isotropic radiation. 

We therefore treat separately radiation emitted (from 
internal heating sources or from external light absorbed by 
the disk) which is mainly isotropic and thermalized, and 
the radiation coming from external sources. The latter can 
also be split into two terms, the attenuated radiation com- 
ing from its external source in a certain known direction, 
and the radiation that has been scattered one or several 
times which should be mainly isotropic. In the first part 
of this section, following Hubeny we neglect the radiation 
by external sources. Later in Sect. 12.41 we will take into 
account correctly the infiuence of external radiations. 

2.3.1. Energy dissipation in each layer 

We first express the temperature as a function of the lo- 
cal production of energy u, such as viscous heating, stellar 
irradiation, disk backwarming, heating by a UV or X ra- 
diation field, or the contribution of the e nergy horizonta l 
transfer. Following the usual notations l)Mihala,s I 
we Hnk u to the source function S and the diffuse radia- 
tion intensity /: 



2.2. Hydrostatic equilibrium — Density distribution 

For a given vertical temperature law, we compute the den- 
sity distribution from the pressure given by the hydro- 
static equilibrium. As the disk is dominated by gravitation 
of the central star, the partial derivative of Gz is a constant 
in the approximation z <ti r, and we use a differentiated 
hydrostatic equilibrium equation for convenience: 



d^P 
dm 2 



1 



dGz 



S{P,T) dz 

The second order equation Eq. is solved with two 
following boundary conditions: 

1. in the disk mid-plane: the first derivative of P is zero 
by symmetry. 

2. in the outer boundary: P tends to a constant low 
value, the pressure in the interstellar medium (here- 
after ISM) . The choice of this limit has almost no in- 
fiuence on the disk structure. 



2.3. Radiative equilibrium — Temperature distribution 

To complete the calculations of the distributions of the 
physical quantities, we have to connect the local tempera- 
ture with the heating sources, remote or local. This is the 
goal of the radiative transfer described in this section. 



471' / Xi/('Ti) [S^{m) — J,j{rn)] dv — p{m)u{rn) (5) 
Jo 

with Jy{m) = j (to, //) d/i (6) 

The conservation of the radiative energy means that in 
an atmosphere layer, the radiated energy, i.e. the emission 
of the layer, minus the radiation coming from the neigh- 
bouring layers Ju{,z), is exactly equal to the dissipated 
energy in the layer: 



1 



J V 1 



(7) 



Integrating over frequencies, one derives from Eq. ®: 



kbB = KjJ + — 
47r 



(8) 



where kb and kj are respectively the Planck absorption 
mean and the usual absorption mean, defined by: 



Kb 



{k^I P)jy J 



(9) 
(10) 



Since -kB = ctbT^, the vertical temperature distribu- 
tion is, by using Eq. JHJ: 



^4 ^ ^ Kj 



J 



47rK ] 



(11) 
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In other words, the temperature in a layer is the sum of 
two contributions: (i) a term proportional to the radiative 
intensity J transferred by adjacent layers and (ii) a local 
dissipation term proportional to u. 

2.3.2. Radiative transfer of the diffuse intensity 

We solve the transfer of the diffuse intensity / to determine 
J. The equation of transfer and its two first moments are: 

d lu 

l^~g^(^' ^) = Xi'iz) [Si,{z) - Iu{z, fi)] (12) 



dz 
dK^ 

dz 



{z) = Xv{z) [Si,{z) ~ J^{z) 
[z) = -Xu{z)H^(z) 



(13) 
(14) 



where J^, Hi, and K^, are respectively the zeroth, first and 
second moments of I^, . 

By integrating over frequencies one gets: 



dH _ u 
dm Att 

-, — - XhH 
am 



(15) 
(16) 



In this equation, the fiux-averaged opacity xh is defined 
by: 

XH^ f {x./p)H,dv/H, (17) 
Jo 

When integrating the first moment equation Eq. Ijl5ll 
over m, one gets: 

H{m) = H{0) - UMe{m) /A-k (18) 

I'M 

with UM = / u(C) dC (19) 



and 9[m) 



1 



UM 



^(C)dC 



(20) 



where U is the total energy dissipation in the upper atmo- 
sphere of the disk per surface unit, and 9 the distribution 
function of the energy dissipation in this atmosphere. 

The symmetry conditions in the equatorial plane for 
the fiux H{m) impHes that H{M) = 0, that is 6'(Af) = 1. 
We therefore get a relation between H{Q) and the energy 
flux dissipated in all atmospheric layers U: 



H{0) = UM/Att 



(21) 



UnHke in a stellar atmosphere, energy is locally produced 
so that the flux is not constant. 

Integrating the second moment equation Eq. Ijlfill over 
TO, one derives: 

K^K{Q) + H{Q) (th-Ath) (22) 
with TH{m) = / XH(C)dC (23) 



(24) 



and ArH(TO) - / XniOmdC 
Jo 



where th and Ath are respectively the flux weighted and 
the "energy dissipation weighted" mean optical depths. In 
a stellar atmosphere, one has K = K{0) + H{0)th- Here, 
Ath represents a correction to the flux-weighted mean 
optical depth due to the decrease of the flux toward the 
mid-plane of the disk. In the large depth approximation, 
we have: dT'^/drn oc H. Since in a stellar atmosphere, 
the flux H is constant, oc Hth, whereas in a disk H 
decreases with th until T"* < H{0)th, hence the corrective 
term Ath- 

We then introduce the Eddington factors. 



f K H{0) 
IK — — and fff = — -— . 
J ■' J(0) 



(25) 



In the large depth approximation, met in thick parts of 
the disk, they are respectively equal to 1/3 and 1/2. 

2.3.3. Formal solution 

Combining Eq. JTll, Eq. and Eq. we get the 

temperature distribution: 



j,4 _ '^J^cff 



K 



TH - Ath 



/k(0) 

fH 



Mkj U 



(26) 



with (JbT^s = UM. 

In the case of strict thermal equilibrium, the intensity- 
weighted and flux-weighted mean opacities kj and kh can 
be substituted by the Planck and Rosseland mean opaci- 
ties KB and KB.- 

Equation Ij26|l means that the vertical distribution of 
the temperature in a disk is similar to its distribution in 
a stellar atmosphere: 



T\t) = 



with two main differences: 



/k(0) 

fH 



(27) 



(i) a local term, proportional to the energy dissipation u 
is added. In a stellar atmosphere no energy dissipation 
occurs in the intermediate layers. The factor propor- 
tional to the density of energy per unit mass u, in 
Eq. (l26|l is also inversely proportional to the optical 
depth of the disk, kjAI. This factor is therefore dom- 
inant in zones where the disk is thin. 

(ii) the optical depth r becomes th — Ath. The term Ath 
is positive but cannot exceed th depending where the 
energy dissipation occurs in the vertical structure. The 
two extreme situations are when the dissipation func- 
tion is a Dirac function located either in the central 
plane of the disk, or at the top of the photosphere. 
In the flrst situation, Ath — everywhere, whereas 
in the second one Ath — th- Above the layers where 
most dissipation occurs, the correction Ath can be ig- 
nored and the disk behaves Hke a stellar atmosphere: 
there is no local energy source and all radiation comes 
from deeper layers. In the central layers, 9{m) sa 1 
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and TH{m) — AT//(m) is constant: most energy dis- 
sipation occurs in upper layers, and, in other terms 
there is almost no radiation coming from deeper lay- 
ers. Therefore the temperature is more or less constant. 
In the particular case where the energy is dissipated 
only in the central plane th — = th and the disk 
photosphere behaves like a stellar atmosphere. 

2.3.4. Several sources of energy dissipation 

In the case where there are several sources of energy dis- 
sipation releasing the energy Uk{m) in the layer m, the 
formal solution given in Eq. ^2<^l can be split as the sum 
of the temperature distributions (see appendix fXl: 



(28) 



with the temperature distribution of source fc. 



^j(^cff)| 



TH - (At//) A; 
Mkj Uk 



MO) 

fH 



(29) 



where the effective temperature of the fc-th source is de- 
fined by: 



aBiTes)t^ / Ufc(C)dC 

Jo 





UkM. 



(30) 



The quantity (AT/f)^ is defined by Eq. lf2ijl where the 
vertical distribution function 9{m) is replaced by the cor- 
responding fc-th energy source. 



UkM 



^^fc(C)dC. 



(31) 



We see that when dealing with several sources of en- 
ergy dissipation, one can treat them separately from the 
analytical point of view (see Sect . 12 .411 . However in the full 
modelization, the field of radiation which defines J(m), 
H(m), fKim) and Jh must be computed globaly. 

2.4. Influence of external radiative sources 



iHubenv I l)l99nl) already took into accour it the effect o f 
the irradiation by the central star (see also lHube nv"l99lV 
However, when addressing the infiuence of the radiation 
from external sources, we can use Eqs. Ij28|l and Ij29|l taking 
into account as an additional source of energy dissipation 
the reprocessing of the incoming radiation. To compute 
the reprocessed energy ua and its influence on {Ath)a, 
we need to solve the radiative transfer of the direct and 
scattered incoming radiation. 

The flux coming from a central source is highly direc- 
tional and cannot be processed as isotropic radiation. For 
convenience, one may split the radiation fleld into several 



components, and consider separate partial transfer equa- 
tions for the individual components. Some of the radiation 
components may be considred as a source of additional 
heating. Only the component containing the thermal emis- 
sion term, k^-B^ can be treated with the formalism devel- 
oped in the previous section. The mean opacities kj, kb, 
Xh, and the mean Eddington factors fx and fn are then 
deflned as appropriate averages over the moments of this 
particular component of the radiation fleld. 

Here, we process separately (i) the radiation from in- 
ternal heating sources or from external light absorbed 
and reprocessed in the disk, a radiation which is mainly 
isotropic and thermalized, and (ii) the radiation coming 
from external sources which is not absorbed. In this latter 
case, we again split the incoming radiation into two terms, 
one attenuated coming directly from its external source in 
a given known direction, and one that has been scattered 
once or several times and that should be mainly isotropic. 

We use the following notations: 

1. Ii,: radiation emitted by the disk itself. (All ther- 
mal flux related quantities are written without super- 
script.) 

2. I^: radiation coming from external sources and scat- 
tered in the disk 

3. /°: attenuated radiation coming from external sources 

4. uy. the energy dissipated locally (e.g. viscosity pro- 
cess) 

5. Ua: the energy coming from the incident photons ab- 
sorbed in the disk and reprocessed 

Th is a pproac h is q u ite similar to Chandrasekhar's 
one JChandrasekhar I Il96flll and was used by 
iMalbet fc Bertout I l|l99l|) . It consists of writing the 
equations of transfer for the total intensity Ii, and then 
deriving equations involving non-local terms /° and If,. 
Since the radiative transfer of Ii, and non-local terms are 
treated separately, the coupling between them can be 
seen either as a loss of energy for the non-local radiation, 
or as a gain of energy for Ii, . The stellar flux can now be 
seen as energy dissipation similar to that due to viscosity. 



2.4.1. Energy reprocessed by the disk 

In the presence of an external source of radiation, the en- 
ergy dissipation term ma Hnked to the reprocessing of the 
radiation coming from external sources and absorbed in 
the disk is directly linked to the speciflc intensity of the 
incoming radiation by: 



p{m)uA{in) = 47r 



+ Ilim, p)) du dfi 



(32) 
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With J° and the frequency-integrated zeroth order 



moments of /° and , 

r+oc ^+1 

J°-y J°di/ with J" = J 



Il{li)dii, (33) 
J"^j^ J^di. with Jl =y^/:(M)d/i, (34) 



and, the intensity weighted opacities, 

p + ca 

/ {Kjp)4dvlj\ 

Jo 

p + ca 

/ {n,/p)Jldv/J\ 
Jo 

UA can be expressed as: 

UA = 4tt{k°jJ° + k'jJ"), 



(35) 
(36) 

(37) 



The energy absorbed ua is therefore proportional to the 
mean incident flux and to the absorption coefficient of the 
medium. 



2.4.2. Radiative transfer of the incoming radiation 

We study the transfer of the radiative intensities /° and 
P upon which our knowledge of ua and {Ath)a de- 
pends. The equations of transfer for the incoming radi- 
ation within the frame of coherent isotropic scattering: 

dl° 

li-^iz,fi) - -x.I°iz,fi), (38) 

dr 

fi-^{z,n) = -Xull{z.,p) +(T^Jl{z,ix) (39) 

where ayJ^{z, fi) accounts for the scattering of the incom- 
ing radiation and auJ^{z,fi) for multiple scattering. 

We define the direction and frequency averaged opac- 
ities, 



X°7 = 



+ 00 



Xh = 



(x./p)J°dz./j"; 


(40) 


OO 


(41) 


OO 

{a,/p)J°d,y/j"; 


(42) 


OO 

{x./p)Kd,^/H'; 


(43) 



where J°, J^, and, H^, HI are respectively the zeroth 
and first moments of /° and II . Therefore the frequency- 
integrated zeroth, first and second moments of and 
namely J", K'^ for direct incoming radiation and J'', 
, for scattered incoming radiation, are linked by the 
following equations: 







dH' 
dm 
dH"" 

dm 



tO 



and 



-k} J" + (7^} J° and 



dm 
dX" 

dm 



x'hH" (44) 
XhH'- (45) 



Therefore, the reprocessed energy, 

,dH° dH^ 



leads to 



UA{m) 

UaM : 
OA{m) : 



-47r 



dm 



dm 



H°{m)+H%m) 
H"{0) + H^{0) 



1 



(46) 

(47) 
(48) 



assuming H^{AI) = and H'^{M) = 0, i.e. assuming that 
the equatorial plane is a plane of symmetry for the source 
of external radiation. This is truly the case for a central 
star or a uniform ambiant medium. 

Using the second moment of the incoming and repro- 
cessed radiative equations, we derive {Ath)a, 



th 



j° [K° - X°(0)] -I- [K"" - X" 



(49) 



where we define respectively and uj^, the mean ratios 
of XH to x% and Xh- 



X//(C)i/°(C)dC=c.°(m) / x^HiOH'iOdC (50) 



XH{OH'iOdC=uJ%m) XH(C)i?^(C)dC. (51) 

If the radiation is absorbed in a geometrically thin layer, 
where the temperature and density are more or less con- 
stant, then « Xh/x'h- ^he ratio of the absorption 
of reprocessed radiation (generally in the thermal infrared) 
to the absorption of incoming radiation (often in the visi- 
ble, or UV, or X-ray domain). In the case of a hot source 
like stellar irradiation or a UV field, this ratio is expected 
to be less than unity. 

Finally the contribution of the incoming radiation 
source to the temperature distribution is: 



TA'(m) 



tj"[i^"(0) 



(^off)! 



fK{0) fx np'^ + K^jJ^ 



(52) 



We show in Appendix El how to compute J^{m), 
J^{m), H'^im), and H^{m) in the case of an external point- 
like source where the irradiation from the central source is 
highly collimated. As demonstrated in Sect. 12. 3. il the tem- 
perature vertical structure of a sum of point-like sources 
is the sum of the temperatures computed for each source 
to the fourth power. 
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2.4.3. Interpretation 

In this section, we discuss in detail the physical interpre- 
tation of the complex Eq. lf52|l . Since we know how to 
add the contributions of various energy sources by using 
Eqs. Ij28|l and Ij29ll . we can restrict ourselves to a point- 
like source located at infinity without losing generality. 
An extended source will then be processed as the sum of 
point-like sources. 

The dependence of the central temperature to the in- 
coming stellar radiation is relatively simple. Since the top- 
most layer emits half upward to the outside and half down- 
ward toward the inner layers, the latter ones receive half of 
the emission. Their effective temperature is then 
Since there is no flux inside (for t ^ 1), the temperature 
is constant and equal to this temperature. This intuitive 
result is compatible with Eq. Ij52ll . As a matter of fact, if 

is the cosine angle of the incoming radiation and since 
K = fi^H, the flrst term in the square brackets is of the 
order of w^/i^. For a slanted incident stellar radiation, we 
have <C 1 and therefore to'^f/' ^ 1. Deep inside the 
disk, the third term, evanescent since the incoming radi- 
ation is absorbed, is also negligible. The brackets reduce 
to /if(0)//ij « 2/3. The central temperature is there- 
fore smaller than the effective temperature Ta by a factor 
^ 2-i/4_ conclusion, half of the absorbed energy is ir- 
radiated by the hot topmost layer and another half by the 
cooler inner layers. 

At the surface, the third term in the square brack- 
ets reduces to 1/{^'^uj^) by using J° — H/fi'^ and is 
dominant. The temperature is therefore much higher 
than the effective t e mper ature. This fact is explained in 
iMalbet fc Bertouti l|l99ll) : the topmost layer, vertically 
optically thin but radially optically thick, absorbs almost 
all stellar light. This layer is su perheated by a factor of 
(l/At°)i/4. IChiang fc Goldreich I fTflflTi) reflned this view 
by multiplying the term by The first factor is 

a pure geometrical effect of the slanted incidence whereas 
the second one comes from the difference between the 
opacities for the visible incoming radiation and the in- 
frared reprocessed radiation. 

2.4.4. Overview of other heating sources 

The fiux radiation from the ISM can be processed on the 
same basis as the stellar irradiation. The main difference 
is that the ISM does not irradiate in a particular direction. 
If we assume that it consists of an isotropic black body ra- 
diation of characteristic temperature Tjsm, we can apply 
Eq. Ij61|l within the two-stream approximation {i.e. with 
H* = l/VS using the definition given in Sect. I3.2"2)l and 
we can demonstrate that the contribution of the ISM radi- 
ation is a constant. In fact, if the ISM is the unique source 
of irradiation, the disk would be isothermal. 

All other radiative sources can be treated exactly as 
the central star. If they are extended, the previous results 
derived for point-like sources are integrated over angles. 



The radial extent of the disk is much larger than the 
vertical one. We can therefore apply the approximations 
of large depth to the horizontal transfer. This hypothesis, 
combined with the LTE approximation, allows us to write 
the horizontal flux as: 



— — 



rp6 



(53) 



37r(pKR) dr 

The local energy gain for the material between r and r+dr 
is 'iTi{Hr{r+dr) — Hr{r)). One then derives the heating per 
unit of mass of material that allow us to self-consistently 
compute the correction to the vertical transfer, in the case 
Hr < H: 

= ^ 54 

p or 

3. Application to T Tauri disks 

The two main heating sources in the case of T Tauri disks 
are the viscous dissipation and the stellar light absorption. 
We apply our general equations in this frame. 



3.1. Approximations 

All the terms t^, kh, etc., are dependent on the structure. 
Moreover fx and /h, are not known a priori. Therefore 
we will seek a self-consistent iterative method able to de- 
termine these quantities. 

We now make two approximations concerning the 
opacities: the diffuse intensity is a Planck distribution at 
the local temperature and the radiation coming from a 
external source is also a Planck distribution at the tem- 
perature of the source T*. We then derive: 

— Kj = Kb, the Planck mean opacity 

— k'j — Kj = K*g, the mean opacity for a black body 
radiation at T* through a medium at temperature T. 

— kh = kr, the Rosseland mean opacity 

Since the scope of this paper is to present analytical solu- 
tions and simple numerical a pplications, we use the gray 
Rosseland opacities given bv iRe]] Li^ l|l994l) . We also 
use the gray approximation = Kg = kr and ignore 
diffusion. Therefore uj^ = that we now call uj 

The fn and Jk values will be close to 1/2 and 1/3. 
These approximations are only valid at large depths, but 
we refine these values in optically thin zones at each iter- 
ation, as well as for uj. 



3.2. Vertical temperature distribution 

3.2.1. Viscous heating at each layer in the disk 

We consider here an energy density per mass unit 
which depends on local conditions. The energy dis- 
sipation per mass un i t due to viscosity is given by 
iFranck. King fc Raine I (jji 



(55) 
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In the case of a thin Keplerian disk we get: 

9 

uv{m) = -r2^j^visc("^) 



(56) 



Since uy(rn) depends vertically only on ^'visc(TO), the 
9v function depends only on the vertical structure of the 
viscosity i^visc and of its mass-average i^visc- 

Up to now there has been no exact description of the 
viscosity induced by the turbulence. In this paper, we have 
considered two different types of viscosity: 

1. The Shakura-Sunyaev prescription stating that the vis- 
cosity is proportional to the local sound speed and local 
height scale: 

i/visc = acsh (X T (57) 

The 9y function is therefore the distribution function 
of the temperature. The term th — {Ath)v, by def- 
inition smaller than th, can be much smaller if the 
temperature peaks at the surface of the photosphere. 

2. The /3-prescription, derived from laborat ory experi- 
ments , and applied to accretion disks by iHure et "all 
l|2nonl) . states that 



(58) 



where /3 is a constant ranging from 10~^ to 10"'^. Here 
the viscosity is uniform along the vertical axis and does 
not depend on the temperature distribution of the disk. 

3.2.2. Heating by the central stellar source 

With the approximations of Sect. 13.11 there is no reason 
to distinguish the two specific intensities /° and I". If we 
call /* + then we f ind an analytical so l ution very 
similar to Eq. (28) of from lMalbet fc Bertoutl Jlflflll) : 



H*{0) 

/a'(0) , fx^yj* 



(59) 



where the quantities Hnked to the star incoming radiation 
are superscripted with *. The only difference is the pres- 
ence of a;* whi c h is n ot a priori equal to 1 as noted by 
iD'Alessio et~an l|l998l) . 

We derive J*(0) and H*{0) from the geometrical prop- 
erties of the system using the approximation described by 
iRuden PolladTI The star is seen by the disk as 

a point-like source. It is valid beyond a few stellar radii; 
in the inner parts this is not the case, nevertheless re- 
pr ocessing is seldom do minant in this region (see Fig. 3 
of lP'Alessio et al.lll998ll . A point-like source equivalent to 
the star must present the same values for J*(0) and H*{0). 
The incidence angle is given by: 



^* = H*{0)/J*{0) 



(60) 



One finally obtains, using the gray optical depth r, the 
same expression (33) of Malbet & Bertout (1991) except 
for the factor ui* already discussed in the previous section, 



MO) 



-r/fj.' 



K- 



(61) 



where the effective temperature of the stellar reprocessing 
is: 

(Teff)i = 4^if*(0)/(7B (62) 

3.3. Equation of state — hydrostatic equilibrium 

The equation of state in the disk is: 



with Cs = 



fcRT 



(63) 
(64) 



First of all, r adiation pressure is ignored in Eq. Ij63|l as 
suggested bv IShakura fc SunvaevTl|l973ll. If the disk is 
assumed to be homogeneous and mixed by the turbulence, 
the mass fraction of gas, /gas, is a known function of T and 
p, and might even be a constant if evaporation is ignored. 
The dust is only a few percent of the mass of the disk, so 
we use /gas = 1 in the present work. 

We also assume that the disk dynamics is dominated 
by the central star, so that it is Keplerian. Then, in the 
thin approximation, Gz is a linear function of z. Combined 
with Eq. Ij63|l . the hydrostatic equilibrium and the local 
height scale become: 



d2p 



2r,2 



dm? 



p 



with Q 



GM* 



(65) 
(66) 

(67) 



3.4. Self-similar accretion model 

The radial distribution of energy dissipation for the ac- 
cretion UY{r) and the surface density S(r) are essential 
parameters for the initialisation step of our code. In or- 
der to be as general as possible, we use a prescription 
with a self-similar ac c retion rate, instead of the model by 
IShakura fc Sunva,ev I l|l973l) which assumes a uniform ac- 
cretion rate over the disk. We assume that some ejection 
process removes material away from the close environment 
of the star within a radius tq. A mor e detailed description 
of this mechanism can be found in iFerreira fc Pelletier I 
(■199,'i) . The accretion rate is parametrized as: 



M(r) 



Moo - 

Vo 

Moc: 



if r < To 
if r >rQ 
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Table 1. Parameters used in the model 





= 2Rq 


M, 


= 0.5 Mq 


T, 


= 4000 K 


TlSM 


= 15K 


ro 


= 100 r. 


Q 


= 0.01 


^min 


= 2.2 Rq 




= 1 


d 


= 450 pc 







where ^ f ^ 1 is the ejection index (see 
iFerreira, & Pelletier and ro the cut-off radius. When 

^ = 0, this model corresponds exactly to the standard ac- 
cretion model bv lShakura fc Sunvaev 

The dissipation per area unit U\{r) and the mass col- 
umn S(r) can then be expressed as: 



= ^— f{r) 



(69) 
(70) 



where /(r) is a factor clo se to unity at large r. Using 
the same method as .Shak ur a fc Sunya.ev I l)l97.j) we ob- 
tain: 



fir) 



l + 2£_ 



(t) 



if r < ro 



2e 



1 



(71) 



3.5. The numerical code 

Our set of equations are coupled and non-linear. They 
need to be solved numerically, within the framefork pre- 
sented below. 

The grid: we use a logarithmic radial grid. For each 
radius we chose a non-uniform vertical grid so that the 
mass column coordinate m{z) is logarithmically sampled. 

Initial conditions: from the radial structure described 
in Sect. 12.11 we compute an isothermal vertical structure. 

Iterative method: Knowing the temperature and mass 
column, the hydrostatic equilibrium is solved for. We de- 
duce the density and pressure distributions. Then the tem- 
perature and optical depth are computed again using the 



development of Sect. 12.3.21 Subsequent quantities are de- 
rived, Hke the kinematic viscosity. The value of M[r) is 
refined so that it matches Eq. ifTflll using the average kine- 
matic viscosity. The mass column grid is computed again 
with the new value of M{r). The latter step is performed 
as many times as necessary to reach convergence (typically 
in ten iterations). 

The parameters used for the computations are re- 
ported in Tabled 

4. Results and discussion 

In this paper, we chose to treat mainly the radiative trans- 
fer in the disk. Even if we have developed the formaHsm 



of the central source heating, we will ignore its contribu- 
tion in a first step and concentrate only on the viscous 
dissipation. In a future paper, we will present the case of 
heating by the central star. Moreover the heating by the 
central star becomes dominant at radii l arger than a few 
astronomical units as shown in Fig. 3 of iD'Alessio et alJ 
l)l998l) so the results presented here will be correct only 
for the inner part of the disk and therefore for the visible 
and infrared thermal domain. 

First, we compute the physical conditions in the disk: 
temperature, density, optical depth. Then we use these 
physical conditions to compute astronomical observables 
like SEDs, millimetric images and interferometric complex 
visibilities and direct images. 

4.1. Physical conditions in the disk 

4.1.1. Temperature and density 

A map of density and temperature for a standard disk 
(^ — 0) with an accretion rate of 10~^ Mq/jv is displayed 
in Fig. The shape of the temperature iso-contours il- 
lustrates that the disk temperature is almost constant in 
the optically thick part, whereas the temperature changes 
quickly with the altitude in the upper layer. As previously 
guessed, the disk appears roughly isothermal over the ver- 
tical dimension in its inner part. 

FigurelSldisplays the variation of the temperature with 
the radius. We verify that the effective temperature, i.e. 
at r = 1, is proportional to the dissipated fiux: it follows 
an r"'^/^ law. We stress that the central temperature is 
much higher than the effective temperature. 

We check that the disk is optically thick even at large 
radii while the geometrical thickness of the disk remains 
much smaller than unity: if we define the disk surface as 
the region of optical depth unity, we find z/r « 0.1-0.2. 

The geometry of the density contours displayed in 
Fig. n] have a different shape depending whether they are 
located in the central optically thick part or in the up- 
per optically thin layers. Figure also displays the disk 
parts which are unstable to convection by computing the 
Schwarzschild criterion. For the time being, the conse- 
quences of this convection are not taken into account in 
our model. 

4.1.2. Column density 

Figure 13 shows the mass column density as a function of 
the radius in the case of a standard a-disk with two values 
of the ejection parameter ^. We take the cut-off ro at 100 
stellar radii, i.e. about 1 AU. Two values of the accretion 
rate have been considered: 10~^ and lO~^M0/yr. As ex- 
pected, the radial structure of S(r) changes in the close 
environment of the star when £,^0. 

Power-laws are often used by observers in the descrip- 
tion of the radial structure of disks, though one can see 
in Fig.01that the column density of an a-disk cannnot be 
described with such a law. However, on restricted domains 




radius (a.u.) 

Fig. 1. Iso-contour maps of a T Tauri disk with M = 
lO~^M0/yr. Dashed lines: lines of equal temperature; 
Dotted lines: lines of equal density; Solid line: limits of 
the optically thick region (r — 1); Gray zone: zone of 
instability with respect to the convection. 




0.1 1.0 10.0 100.0 

radius (a.u.) 

Fig. 3. Column density in an a-disk for two different ac- 
cretion rates. Solid lines: ^ = 0; Dashed lines: S. = 0.5 



radius (a.u.) 

Fig. 2. Radial distribution of the temperature in the cen- 
tral layer and of the effective temperature (located at 
r = l). 




0.1 1.0 10.0 100.0 



radius (a.u.) 

Fig. 4. Column density in a disk with different viscosity 
prescriptions. Solid lines: a = 10~^, 10~^; Dotted lines: 

(3 = lo-^ 10-3 



of radius, the mass column in an a-disk can be described 
by the equation 

S(r) « AU)-'' kg.m'^ for r, < r < n+i (72) 



where the coefficient qi ranges from -0.4 to 1.5 and Sj is 
a constant. These coefficients are given in Table for two 
values of the accretion rate. 

We show in Fig. the influence of a on the struc- 
ture of the disk. These results are compared with the ones 
obtained with the /? viscosity prescription. The general 
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Table 2. Parameters computed from the model to ap- 
proximate the surface density by a power law by pieces, 
sorted by decreasing order of radius domain. 



10+4 - 



M = lO-'^Mo/yr 
Si qi Ti 


M = lO-*M0/yr 

Qi 


1 10" 1.5 40 
410^ 0.0 9 
1.510* 0.6 1 
1.510* 1.2 0.11 
410=5 -0.3 0.06 


2 10* 0.0 67 
2.510=5 0.6 5 

910=5 1.4 0.75 
1.510*5 -0.5 0.45 



behaviour is somewhat different except for the 1 — 10 AU 
region. Therefore observing inside 1 AU or outside 10 AU 
could be decisive in determining the viscosity law. 

4.2. Astronomical observables 

One of the aims of this work is to compute physical con- 
ditions around a young star so that we can calculate the 
electromagnetic field emerging from this region. This emis- 
sion can then be analysed by several types of instruments: 
photometers, spectrographs, imaging cameras, interferom- 
eters, etc. The next sections give examples of astronomical 
observables that can be used to constrain our knowledge 
of the circumstellar environment. 

4.2.1. Spectral energy distribution 

The spectrophotometry technique gives access to the spec- 
tral energy distribution. We integrate along the fine of 
sight the contribution of each disk layer to the emergent 
fiux. In the examples taken for Fig. the SED is sig- 
nificantly altered in the range 1-5 fim with the self-similar 
accretion model. The inner hot regions being depleted, the 
short wavelengths are less present. We assumed an overall 
extinction of Ay = 1 in the visible and an inner radius 
of the disk of rmin — 1-1 stellar radius. However, Fig. |3 
shows that the SED is significatively different from the 
one expected from a standard model with a large inner 
gap. 

Figure shows the influence of disk flaring on the in- 
termediate wavelength range of the SED. The difference 
between the two curves representing the SED of a 45°- 
inclined disk lies in the location of the photosphere, i.e. 
the location of the disk layer that emits most of the emer- 
gent flux. In the flaring case, the flux comes from the layer 
where t — 1, whereas in the so-called no flaring case, the 
flux comes from the disk equatorial plane but with the 
same temperature as the effective temperature. This lat- 
ter case corresponds to the standard radial model where 
the vertical structure is not taken into account. Even with- 
out the reprocessing of the stellar photons, the effect of the 
vertical extension of the disk atmosphere is therefore not 
negligible at intermediate wavelengths. 



T3 

N 

E 
o 



> 
w 

E 





100.0 



radius (a.u.) 



Fig. 7. Emergent millimetric flux from a pole-on disk 
for several accretion rates and viscosity prescriptions. 
Upper curves: M — lO~'*M0/yr; Lower curves: M — 
lO~®M0/yr; Solid lines: a = 10~^; Dashed lines: (3 = 
10-4. 

4.2.2. Millimetric innages 

Radio observations can help us discriminate between the 
two viscosity prescriptions, and especially radio observa- 
tions in the outer disk, since they are sensitive to the mass 
column. The disk being optically thin in the millimetric 
continuum, the flux is given by: 



-{-OQ 



T{z)p{z)dzK ST 



(73) 



Figured compares the millimetric emissivity for an a-disk 
and a /3-disk with the same characteristics. At interme- 
diate radii, the two emissivities are more or less equiv- 
alent, however at larger radii the emissivities can follow 
somewhat different power laws. is proportionnal to 
7.-3/2 £qj. g^j^ a-disk over the full range, whereas a /3-disk 
presents a variable power-law exponent. In the domain 
where the disk is thermalized with the ISM, i.e. at large 
radii, Fu oc r~^/^. 

4.2.3. Interferometric observables 

Probing the inner parts of disks requires high angu- 
lar resolution techniques that can be achieved with 
presen t or soon-to-be operated interferometric instru- 
ments l|Ma,lbet fc Bertoutlll99,'j ) . The more massive disks 
(M > lO"'^ Me/yr) should be resolved by the VLTI and 
its lOOm-baseline, but we should also detect non-zero clo- 
sure phase. 

Figure |H1 shows the visibility of a disk as a func- 
tion of the projected baseline. Since the image is not 
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10-10 




100.0 




wavelength (10"^m) 



Fig. 5. Spectral energy distribution {XF\ vs. A) for differ- 
ent disk models Solid lines: standard disk; Dashed lines: 
disk with ^ = 0.5; Dashed & doted lines: standard disk 
with an inner hole of 8r* . 



100.0 



wavelength (1 0"^m) 



Fig. 6. Spectral energy distribution {XF\ vs. A) of the 
disk model as a function of inclination. Solid lines: pole- 
on disk. Dashed lines: inclination of « = 45°. Dotted li- 
ness: same inclination without the flaring effect (see text 
for details). 



centro-symmetric, the visibility depends on the position 
angle. Figure El displays the closure phase as a function 
of the position angle for different telescope combinations 
in the VLTI for a disk with an important accretion rate 
(M = 10~^ Mq/yt). As a comparison, the closure phase 
is displayed for a disk with moderate accretion rate with 
the widest combination of the VLTI auxiliary telescopes. 
Because of the natural flaring, a disk that is not seen pole- 
on does not appear centro-symmetric (as can be noticed 
in Fig. nTijl . therefore leading to a non-zero closure phase. 

4.2.4. Direct images 

Synthetic infrared images of the disk thermal emission are 
displayed in Fig. ^1 with a fleld of 50 mas for a distance 
of d = 450 pc. As explained above, they are not centro- 
symmetric because of the flaring. So far, we have only ac- 
cess to interferometric observations at this scale, so they 
are only relevant as a step to produce interferometric ob- 
servables. 

5. Conclusion 

We have presented a new method to model the vertical 
structure of accretion disks. This method leads to analyt- 
ical formulae for the temperature distribution which help 
to understand the behaviour of the radiation propagated 
inside the disk. Our model includes two sources of energy 
dissipation: viscous heating and reprocessing of external 
radiation like that emitted by the central star. We have 
applied these analytical results to the case of T Tauri disks 



in a variety of conditions showing that the method is ver- 
satile: 

— with a modifled Shakura-Sunyaev prescription 
for a disk with a non -constant accretion rate 
jFerreira, fc Pelletier lITflfl^ 

- with a viscosity prescription follo wing either the st an- 
dard a or the /3 one proposed bv lHure et alJ 1)200(1) 

We are able to simulate the conditions of temperature 
and density in any part of the circumstellar environment 
and to compute astronomical observables like SED, op- 
tical and millimeter images or interferometric visibilities. 
This code could be used also as starting conditions in a 
Monte-Carlo multiple scattering code in order to derive 
polarization maps. 

When completed, this code will offer to observers a tool 
based on physical parameters to interpret their observa- 
tions better than ad-hoc models based on power laws often 
used today. 
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Appendix A: Vertical structure of the 

temperature with several energy sources 

Each source of energy locally contributes to the heating Uk 
(see Sect. l2.3.i)l . The total local heating is then u = J2k 
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Fig. 8. Visibilities from a disk model as a function of 
the baseline. Gray: M lO^'^Mo/yr; Black: M = 
lO~^M0/yr. For a given baseline the visibility varies 
with the position angle; the domain between the mini- 
mum and maximum values is filled. 
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Fig. 9. Closure phases as a function of the position 
angle for a triplet combination of the VLTI auxiliary 
telescopes. Solid lines: M = lO~'*M0/yr. Dotted line: 
M ~ 10~^ Mq/jt with the largest baseline. 




and the energy dissipated per unit of disk surface is then 
U = Uk with 



UkM = 



M 



UkiO dC. 



(A.l) 



We now write the vertical distribution of energy dissipa- 
tion, 6* as a function of the different vertical distributions 
of the energy source, 9k: 



e{m) 



with 6k{m) 



u{C) dC /{UM) 



UkiOdC/ (UkM) 



(A.2) 



Similarly for Ath: 

ATnim) 

with {ATH)kim) 
Since, 

u{m) 
U ^ 



^(C)rH(C)dC 



(A.3) 



, U Uk 



(A.4) 



and J2ki^k/U) = 1, we can rewrite Eq. Ij26|l by replacing 
each term in the brackets by a sum over the k sources: 



rp4 



E 



Uk 



4:^3 Jk V ^ 



TH - {ATH)k 
!k Uk 

' Mkj Uk 



(A.5) 
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If we define the effective temperature of the the k-th 
energy source, 

cTB{T,s)t ^ UkM, (A.6) 
then Eq. lf2fijl can be written as the following sum: 



rp4 



E 



^./(^cff)fc 



th - (Ar//)fc + 



Ih 



Mkj Uk 



(A.7) 



Therefore the vertical distribution of the temperature for 
a disk with multiple energy dissipation sources is the same 
as the sum of the individual temperature distributions us- 
ing the global radiation in the disk, i.e. with J, H , K, th 
computed from the global temperature distribution. 

Appendix B: Incoming radiation field 

We show in this appendix how to compute J°(m), J^{m), 
H^{m), and H^{m) in the case of an external point-like 
source. The intensity of the source is J", giving us a 
boundary condition to solve the equations of transfer lilill 
and If we define the optical depths, 



we obtain: 



Xj(C)dC 

/o 
Jo 



J°{m) ^ J(0)e-^''(")/'''' 



(B.l) 
(B.2) 

(B.3) 
(B.4) 



where A is a commonly used operator (cf. iMihala.s lll97Stl 
defined by: 

Ar{f)^\£^ E,\t-T\f{t)dt (B.5) 

/ + 00 
u-"e-"*du. (B.6) 

The terms H^{m) and H^{ra) can be computed in the 
same way as J'^{m) and J^{m): 



W{m) = $,s(„,) ( ^ J° 



(B.7) 
(B.8) 



with 



$,(/) = 2/ E2\t~T\f{t)sgn{t-T)dt (B.9) 
Jo 

where sgn(x) is the sign function. 
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